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ABSTRACT 

Our current knowledge of the rotation curve of the Milky Way is remarkably 
poor compared to other galaxies, limited by the combined effects of extinction 
and the lack of large samples of stars with good distance estimates and proper 
motions. Near-future surveys promise a dramatic improvement in the number 
and precision of astrometric, photometric, and spectroscopic measurements of 
stars in the Milky Way's disk. We examine the impact of such surveys on our 
understanding of the Galaxy by "observing" particle realizations of nonaxisym- 
metric disk distributions orbiting in an axisymmetric halo with appropriate errors 
and then attempting to recover the underlying potential using a Markov Chain 
Monte Carlo (MCMC) approach. We demonstrate that the azimuthally averaged 
gravitational force field in the Galactic plane — and hence, to a lesser extent, the 
Galactic mass distribution — can be tightly constrained over a large range of radii 
using a variety of types of surveys so long as the error distribution of the mea- 
surements of the parallax, proper motion, and radial velocity are well understood 
and the disk is surveyed globally. One advantage of our method is that the tar- 
get stars can be selected nonrandomly in real or apparent-magnitude space to 
ensure just such a global sample without biasing the results. Assuming that we 
can always measure the line-of-sight velocity of a star with at least 1 kms~^ 
precision, we demonstrate that the force field can be determined to better than 
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~1% for Galactocentric radii in the range i? = 4 — 20 kpc using either: (1) small 
samples (a few hundred stars) with very accurate trigonometric parallaxes and 
good proper-motion measurements ( uncertainties Sp^tn ^ 10 /xas and 5^ < 100 
/xasyr"^ respectively); (2) modest samples (~ 1000 stars) with good indirect 
parallax estimates (e.g., uncertainty in photometric parallax (5p,phot ~ 10%-20%) 
and good proper-motion measurements (5^ ~ 100 /iasyr"^); or (3) large samples 
(~ 10^ stars) with good indirect parallax estimates and lower accuracy proper- 
motion measurements {5^ ~ 1 masyr"^). We conclude that near-future surveys, 
like SIM Lite, Gaia, and VERA, will provide the first precise mapping of the 
gravitational force field in the region of the Galactic disk. 

Subject headings: dark matter - Galaxy:disk - Galaxy:fundamental parameters 
- Galaxy: kinematics and dynamics - methods: data analysis - surveys 



1. Introduction 



Observations of the motions of stars and gas in galaxies tell us that they contain many 
times more mass in encompassing dark matter halos than in their stellar components (e.g., 
KentI 119871 ). However, exactly how this dark matter is actually distributed in galaxies is 
still of some debate. For example, while simulations of cold dark matter halos forming in 
an expanding universe seem to ge nerally converge on adensity distribution that can be 
represented by a universal formula (INavarro et al.lll996l . 119971 ). the shape an d radial profile 



of the inner parts of dark matter halos are still uncertain (see discussion in iNavarro et al. 
20041 : iHavashi et al.l[20ol . 



Of course, baryons are expected to complicate the elegant simplicity of the picture of 
dark matter halos painted by pure A^-body simulations. Gas radiates away energy to sink 
toward the centers of the dark matter halos where it can contribute significantly to the 



gravitational potential. This process can cause the 
further in response (as reviewed i n 



spher ical in shape (jPubinskil 1 19941 : iKazantzidis et al 



Gnedin et al. 



background dark matter halo to contract 



20041) and evolve from triaxial to more 



2004; Bailin et al. 2005: Abadi et al. 



20091 ). On the other hand, stellar bars at the centers of galaxi es can transfer a ngular momenta 
to their host halos, flattening their central density cusps (ISellwoodI l2006l ). The decay of 
satellite galaxies and substructure can also flatten the central density cusps. 

Ultimately, we want to be able to distinguish between dark and luminous contributions 
to the distribution of matter throughout galaxies. Stellar disks provide some of the cleanest 
probes of matter distributions, with stars moving on near circular orbits. Nevertheless, there 
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remains the tricky problem of decomposing a disk galaxy potential into disk, bulge and halo 
contributions in order to isolate the form of the dark matter distribution. One approach to 
this dilemma has been to look at low-surface-brightness galaxies, which are expected to be 
dominated by dark matter, yet even in these case s the results have been controversial and 
ambiguous (see lde Blokll2005uHayashi et al.ll2004l . for two opposing views) 



It is striking that our own Milky Way galaxy has as yet contributed little to these 
debates. After all, this is the one galaxy we can expect to study star-by-star with very 
high resolution in three-, four- or even six-dimensional phase space. So far, three effects 
have hampered these ambitions: first, our lack of accurate distance measurements to stars; 
second, our lack of accurate proper motions of stars; and third, our inability to see across 
the Galactic disk because of dust absorption. Because of these we do not yet have the solar 
circular speed to better than 10%, the disk scale length to better than 20%, or an accurate 
asses sment of our own Galaxy's rotation curve beyond the solar circle (see lOUing fc Merrifield 
19981). We have only fairly recently bec ome convinced of the barred nature of the Milky Way 
tehtz fc Spergell Il99ll: IWeinberel I1992I ) and are unsure whether we live in a fiocculent or 
grand design spiral (lQuillenll2002l ). 



The Hipparcos Space Astrometry Mission revolutionized our understanding of the solar 
neighborhood by comp iling 1 milliarcsec level astrome try of 120,000 stars. Using this data 
Creze et al.l (119981 ) and lHolmberg fc FlynnI (120001 . l2004l ) measured the local matter density in 
the disk (the Oort limit) to be ~ O.IMqPc"^, a value that leaves little room for any significant 
contribution from disk dark matter. Such an explicit decomposition of baryonic and dark 
matter contributions to a disk potential is impossible in external galaxies. iFlynn et al.l (120061 ) 
used these results to estimate the local surface mass-to- light ratios {M/L) for the Galactic 
disk of {M/L)y = 1.5 ± 0.2 MqLq^ and inferred that the Milky Way is under-luminous 
by about la with respect to the Tully-Fisher relation; if the rotation speed announced by 
Reid et al.l (120091 ) is correct this discrepancy is even more significant. While these studies 
demonstrate the importance of large-scale, systematic Galactic studies to understanding 
galaxies in more detail, Hipparcos' distance horizon was about 100 pc (distances of 10% 
accuracy) so it could not map the distribution of the mass in the Galaxy beyond the solar 
neighborhood. 

Three innovations in observations promise to dramatically improve our understand- 
ing of the phase-space structure of our Galactic disk: (1) large-scale photometric surveys, 
both existing (the Two Micron All Sky Survey (2MASS) and the Sloan Digital Sky Survey) 
and planned (PanSTARRS and LSST), toget her with methods of deriving accurate pho- 
tometric parallaxes for stars in these surveys ( Majewski et al. 2003; Juric et aD bood ) . (2) 
high-precision (few to lO's of /xas) astrometry from radio observations of masers e.g., VERA 
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(iHonma et al.ll2000h. V LBA flReidll2008l : iHachisuka et al.ll2009h and the European VLBI Net- 
work (IRygl et al.ll2008l ) and optical observations of stars (NASA's SIM Lite — Space Interfer- 
om etry Mission Lite and ESA's GAIA — ^ Global Astrometric Interferometer for Astrophysics, 



see 



Unwin et all booll : iPerrymanI |2002| ) : and (3) large-scale, hi gh-resolution spectros copic 



surveys, such as the ongoing Radial Velocity Experir nent (RAVE: ISteinmetz et al.ll2006l ) and 
the SEGUE project of the Sloan Digital Sky Survey (IBeers et al.ll2004j ) as well as the planned 



Apac he Point Observatory Galactic Evolution Experiment (APOGEE; lAUende Prieto et al. 
20081 ) . HERMES instrument for the Anglo Australian Telescope and Wide Field Multi 



Object Spectrograph (WFMOS) for the Gemini telescope. It is clear that any or (better 
yet) all of these advances will significantly improve our knowledge of the Galaxy. What is 
unclear is the relative contribution of each type of survey: how uncertain and/or biased will 
our mass estimates be if one (or more) dimensions of phase space remain unmeasured? How 
far across the Galactic disk do we need to probe in order to construct its rotation curve 
confidently? To what extent can measurement errors be compensated for by using large 
numbers of stars? Our study represents a first step toward addressing these questions. 

Here we describe a general method to recover the underlying potential of the Galaxy 
from photometric, astrometric, and spectroscopic surveys of disk stars (Section [2.2l) . We test 
our method by constructing nonaxisymmetric particle disks orbiting in a given potential, 
simulating observations of these particles with varying degrees of accuracy, sample size, 
and disk coverage (Section 12. ip and examining how well the underlying potential can be 
measured. We present the results of applying the recovery routine to our "observed" data 
sets in Section [31 discuss the implications of these results for future surveys in Section HJ and 
summarize our conclusions in Section [5l 



2. Methods 

2.1. Particle Disk Realizations 

This paper focuses on how astrometric determinations of the motions of disk stars in 
the Galaxy can best be utilized to measure the total potential in which they are moving: we 
neither attempt to disentangle the disk and halo contributions to the potential nor model 
motions perpendicular to the Galactic plane, although the methods that we describe in this 
paper can easily be extended to these tasks. Moreover, we do not address the use of radial 
velocity-only surveys. We explore the power of astrometric measurements to measure the 



This actually presents about SIM PlcLuetQuest instead of SIM Lite 
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Galactic potential by using an approximate, parameterized kinematical model to generate a 
stellar sample, adding observational errors and then attempting to recover the parameters 
of the model. 

The positions and motions of particles in our models are generated from analytical for- 
mulae derived for axisymmetric disks perturbed by spiral arms. We stress that our use of 
approximate analytical formula (e.g., from epicycle theory) does not compromise the valid- 
ity of our results so long as the same approximate formulae used to generate the particle 
realization are used to recover the parameters of the model from observations of the real- 
ization. In other words, the tests described in this paper provide an accurate assessment 
of the validity of our method so long as accurate physics is used when analyzing the real 
data. Of course, a potential problem for applying this method — or any method based on 
parameterized models — is that the results may be misleading if the parameterized models 
are not an accurate description of the real Galaxy. 

The total potential is written as a sum of axisymmetric $ and spiral arm $s terms 



where R is the Galactocentric radius and is the azimuthal angle in the disk, measured from 
the Sun-Galactic center (GC) line and increasing in the same direction as Galactic rotation. 

The particles are assumed to be drawn from an underlying axisymmetric distribution of 
number density Ssym(/2), whose response to the spiral arm potential perturbation T,s{R, 0, t) 
is calculated in the linear regime, to give a total number density: 



The motions of the particles in the underlying potential are chosen to maintain the 
number density distribution: the mean radial (U/j) and azimuthal {v^p) speeds are given by 



where (w_R,sym, v<j>,sym) a^c the mean radial and azimuthal speeds, respectively, set by the 
gravitational potential of the axisymmetric disk, and (f/js, v^jis) are additional perturbations 
to the mean due to spiral structure (see Equations [T4lfT9l below) . 

The formulae, adopted functional form and parameters for axisymmetric and spiral arm 
terms in Equations dH) - (HI) are described in Sections 12.1.11 and 12.1.21 respectively. The 
dashed lines in Figure [1] show the results for our standard model (hereafter, the INPUT 
model), and the dots show the velocities for a sample of particles in the range (p = ±7r/16 
radians generated from this model without errors. 



$total = $(i?) + $s(i?,0,t) 




S = Ssy^(/?) + S,(i?,0,t). 



(2) 



VR = VR^sym{R) +VRs{R,(f),t), 
V4, = V^,sym{R) +V^s{R,(l),t), 



(3) 
(4) 
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2.1.1. Axisymmetric potential, number density and motions 



The realized disks are zero-thickness and exponential in Galactocentric radius R: 



Ssym(-R) = Soexp - 



h 



(5) 



where Eq is the central value and h = 3.0 kpc is the scale length of the disk number density. 

We work in terms of a spherical mass distribution for simplicity, even though the actual 
mass distribution is certainly flattened, because we are only modeling the potential in the 
disk midplane. The combined disk and halo mass distribution is rep resented by a He rnquist 
function (of total mass Mh = 10^2 and scale length a = 22 kpc, lHernquistlll990r ). 



M{R) 



MhR^ 



{R + af 

where M[R) is the mass enclosed within radius R. The corresponding potential is 



$(i?) 



(R + a)' 



(6) 



(7) 



The circular velocity in this potential is calculated via 



t^drc(-R) 



R- 



GM ^JGMhR 



dR \ R (R + a) 
The radial velocity dispersion is assumed to follow 

Rf] — R 



(9) 



where a^^Q = 25 kms~^ is the radial velocity dispersion at the Sun (taken to be at i?o = 8 kpc 
from the GC). Our potential recovery algorithm determines independently of h, i.e., it does 
not assume that = h. If the disk is self-gravitating, the shape of the velocity eUipsoid is 
independent of radius, and the disk thickness is ind ependent of radi us (as observed in external 
galaxies) then we expect h„ = h (see details in iHernquistI Il993l ). Thus our input model 
assumes h„ = h = 3.0 kpc. It i s not certain that this assumption is tr ue for the Galactic 
disk, since estimates for h (e.g., lOiha. Bienavme. Robin &: Mohanlll994l . found h = 2 — 2.5 



Lewis fc Freeman 



19891 . found K = 4.37 ± 0.32 kpc 



kpc for the old disk stars) and (e.g., 
in a study of old disk K-giants) are not same. 

The azimuthal velocity dispersions a,^ are assigned according to epicycle theory. 



a^(R) = (Tr 



(10) 
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where Q = Vdrc/R is the angular velocity of a circular orbit at radius R and k is the epicyclic 
frequency given by 



2,„, 3t;2.^, GM G dM 



The mean azimuthal motion f</,,sym due to the axisymmetric potential is given by the asym- 
metric drift equation, 

where the second equality was derived using Equations and (Q. The mean radial motion 
VR,sym is zero. 



2.1.2. Spiral arms 



The m-armed spiral potential is of the form 



<^s{R,(p,t) =<i>a{R) COS 



m{(j) — Qpt) + clog 



R 



8 kpc 



+ 



(13) 



where $a is the amplitude of the spiral potential, Qp is its pattern speed and k = c/R 
is the radial wave-number which is related to the pitch angle, 6, by cot^^ = c/m. For 
the remainder of the discussion we consider the case of an m = 2 spiral, of amplitude 
^a2{R) = 200(kms~^)2, rip = 1.4 kms"^ kpc-\ pitch angle 9 = 15° and phase (po = 105°. 
This potential perturb ation would be generated by spiral arm of mass density Q.9MqPc~^ 
(from Equation 6.30 of iBinney fc Tremaind 120081 ) or about 10% of the observed local disk 
surface density. 

Note that these parameters are specifically chosen so that any resonances lie outside 
our survey region; in particular the pattern speed is much lower than typical (very uncer - 
tain) estimates of the pattern speed in the Galaxy (jPebattista. Gerhard fc Sevensterll2002l ). 
Moreover, with this pattern speed the entire observable disk lies inside the inner Lindblad 
resonance, a region in which self-consistent spiral waves normally do not propagate. This 
oversimplification will need to be addressed in future work. 

The linear-theory prediction for the velocity response of t he system to the irnpose d 
potential is given in the tight- winding or WKB approximation by (IBinney fc Tremaindl2008l ) : 



VRs{R,(p,t) 
V^s{R,(l),t) 



VRa{R) cos[m( 



R 

' Vtpt) + clog — + 0o], 
R 

v^a{R) sin[m(0 - Vtpt) + clog — + 0o], 



(14) 
(15) 
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where 



VRa{R) 



m{Q — Qp) 

A 
2B 

i;$ P 



and 



B 



djnR) 
dR 



2 dR 



(16) 
(17) 
(18) 

(19) 



In the above equations F is the reduction factor given by Equation (6-63) in lBinney &: Tremaine 
(120081 ) and the function B{R) is equal to Oort's B constant at R = Rq. 



The number density perturbation due to the spiral potential of Equation (fT3l) is given 



by 



T,s{R,(p,t) = Sa(-R) cos 



m(0 — Qpt) + clog 



R 



8 kpc 



+ 



(20) 



where the amplitude of t he response number density is given by Equation (6-59) in 
Binney fc Tremaind (120081 ): 



m{Q - Qp)T.a + kT,oVRa = 0. 



(21) 



2.1.3. Rendering and "Observing" the Particle Disks 

The stellar sample generated from our particle disks simulates the collection of data by 
a targeted astrometric study, such as might be performed with SIM Lite. In particular, the 
target stars are selected based on estimated photometric parallaxes, then the trigonometric 
parallaxes of these targets are observed with appropriate errors. We anticipate that our 
results will be broadly applicable to global astrometric surveys (e.g., the GAIA mission), as 
well as smaller samples currently available from ground-based surveys (e.g., VERA). The 
error distribution function would differ in these cases, but we defer more detailed discussion 
of these alternative applications to Section HI 

We use the following steps to generate our simulated data. 

1. The stars in the disk are generated over a radius range from i? = kpc to 25 kpc, dis- 
tributed according to the number density distribution of Equation ([2]). These stars are 
viewed from a point moving (in Galactic rest-frame coordinates) with radial and az- 
imuthal speeds (wi?,view, ^^0,view) = (^^r.lsr, ^^^.lsr) = {vr,v^) at i?o = 8 kpc from the GC 
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(i.e., at the Sun) and each star is assigned a Galactic longitude / and parallax p. Here 
we assume v^iew = vlsr where vlsr is the velocity of the Local Sta ndard of Rest, be- 



cause the peculiar velocity of the Sun relative to the LSR is known (jPehnen &: Binney 



19981 ). The proper motion n and hne-of-sight velocity vios of the stars are assigned 
based on (fi?,^,/,) drawn from Gaussians with means (yR,V(f,) given by Equations ([3]) 
and (jlj) and dispersions (ci/j, cr^) given by Equations ([9]) and (fTOj) at each position. 

2. The stars' parallaxes are scattered about their true values to give observed photometric 
parallaxes Pphot by drawing from a Gaussian distribution of mean p and dispersion 
5p,phot- This observational uncertainty is assumed to be 15% (i.e., 5p,phot = 0.15p) in 
most of the following work except in Section (see also details in Section [XT]) . This 
step is intended to mimic the effect of drawing a sample from a set of standard candles, 
with absolute magnitudes in a well-characterized range. Errors in I are negligible. 

3. In our recovery algorithm, we parameterize the potential using characteristic masses 
Mj = M{Ri) defined at eight discrete radii Ri uniformly spaced from i?i = 4 kpc to 
Rs = 20 kpc (see Section [2.2.31 for details). In order to have adequate constraints 
on each Mj, sample stars are selected based on observed Pphot (i-e., including the 
uncertainty introduced by the distribution of absolute magnitudes of the standard 
candles) so that there are equal numbers of stars in the seven bins between the Ri. 
We impose the following additional restrictions on the sample, while retaining the 
constraint that there should be equal numbers of stars in each bin. (1) Sample stars are 
selected only between Galactocentric radii of Ri and -Rmax = min(18kpc, Rs — O.lTrf), 
where d is the distance from the Sun. The d dependence in -Rmax is introduced to 
minimize the number of stars that scatter into our sample from outside our intended 
survey region since this could result in a systematic bias in our estimate for Mg. The 
dependence is tuned to the scale of the uncertainties: 0.17d corresponds to the la 
uncertainty in distance due to a 15% error in photometric parallax. (2) No sample 
stars have d > 20 kpc, to avoid stars with large distance errors. (3) Stars behind the 
Galactocentric circle at -Ri (i.e., / < /min = 30° and d > 8 kpc cos /mm) are also excluded 
because of the typically strong extinction. The observed spatial distribution (i.e., using 
Pphot to find location in the disk) of a sample of stars selected in this manner from the 
exponential disk is shown as black dots in Figure [H The true spatial distribution of 
this sample of stars is also shown as gray filled circles. Note that because these target 
stars are selected based on Pphot rather than p, some fraction of our sample actually 
lies outside our intended survey region. 

4. The observed trigonometric parallax, proper motion and line-of-sight velocity {ptri, A^o? 
■^los.o) for our sample are assigned by drawing from Gaussian distributions of mean 
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(p, /i, fios), given in step 1 above, and dispersion ((5p,tri, (^/i, <5ios) where Sp^tri and S\os are 
constant and 5^ depends on 5p tri (see details in Section [311) . This step is intended to 
mimic a targeted astrometric mission with ground-based spectroscopic follow-up, ob- 
serving the sample stars with integration times tailored to achieve a constant accuracy. 
The faintest (i.e., most distant and smallest p) stars in any sample may have more 
accurate photometric parallaxes than trigonometric parallaxes. 



2.2. Recovering the Underlying Potential 

The recovery program uses a Markov Chain Monte Carlo app roach (MCMC, as de- 

and summarized below 



scribed in lGilks. Richardson fc Spiegelhalteiill996l : I Verde et al.l 12003 



in Section 12.2.11) to find the maximum and the shape of the likelihood function 



N 



L{^) = n ^(Ptri, f^l <s,obthot, X), (22) 



for a sample of size as the parameters of the model x are varied. Here P{pl^i, fil,vl^^^\ 
Pphot'^*)^) conditional probability (derived in Section 12.2.21) of observing a star to 

have trigonometric parallax, proper motion and line-of-sight velocity (ptri, yUo, ^^los.o), given 
its observed photometric parallax Pphot and Galactic longitude /, and underlying Galaxy 
model parameters x (see Section [2.2.31) . 

P represents the conditional probability of observing a star's kinematical properties at 
a particular position in the disk. Hence, although we need to have an appropriate model 
of the intrinsic stellar spatial distribution from which we are selecting the sample (in order 
to understand the likelihood of finding a star of given Pphot and ptri in the initial random 
survey), we have complete freedom in specifying how we select the stars in our targeted 
sample. This means we can choose to distribute our tracers to regions of the disk that we 
are most interested in resolving. In our case (as noted in Section [2.1.31) . we select equal 
numbers of stars with observed Galactocentric radii (based on the photometric parallax) in 
each of the seven bins between the Ri, rather than simply taking a random sample, and this 
allows us to explore the outer disk in greater detail. 



2.2.1. The Markov Chain method 



In a single step of a MCMC run, the likelihood I/j,prop (Equation [22]) is evaluated for 
the model parameters Xj^prop proposed at step i in a chain and compared with from the 
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previous step. If Li,prop > qLi^i, for a random number q between and 1, the proposed 
parameters will be adopted for this step (xj = x-j prop). Otherwise the parameters from 
the previous step (xj = Xj_i) are kept. Proposed parameters for the next step (xj+i) are 
generated by adding a vector of small chan ges to x^. These ste ps are accumulated until they 



satisfy the convergence criteria outlined in IVerde et al.l (120031 ) . 



The beauty of the MCMC method is that the distribution of the accepted steps follows 
the shape of the likelihood function in parameter space. This property of the method means 
that the chains of steps themselves can be exploited in two ways: first they can be used to 
derive the optimal directions and sizes of steps for exploring parameter space (i.e., to get 
better acceptance ratios and faster convergence); and second, they can be used to determine 
best values and confidence intervals for each parameter. 

In this work, a "test" MCMC is first run with step sizes for parameter changes es- 
timated from simple intuition. The covariance matrix of these preliminary chains is then 
constructed, and the eigenvalues and eigenvectors of the matrix are used to estimate optimal 
size and direction of parameter vectors of the steps made in the following actual MCMC runs. 
This refinement is particularly important when the parameters have strong correlations (see 
Section [3]) . 

The "best" values of parameters presented in all figures and tables are taken to be 
the mean x from the chains, weighted by the likelihood. The 1-a error bars represent 68% 
confidence intervals. 



2.2.2. Estimating the likelihood of a given parameter set 

The likelihood in Equation fl22|) is the product of factors -P(ptri, A^o, ^^ios,o|Pphot, ^ x) for 
each star in the sample, which is the conditional probability of making observations of 
trigonometric parallax ptri, proper motion = /Xq and line-of-sight velocity v\os = fios.o given 
an observed photometric parallax Pphot, along Galactic longitude /, 

P{n ■ U V. \V. / Xl - :^mot,Ptn,/io,t'los,o|/,x) 
lAri) Po; t^los,o|Fphotj XJ — 



P(Pphot|^,x) 

-^(PphotjPtri; /^o, "i^los.ol^, x) ^^^^ 



/o°° dptri jToo (^l^o JZo dVlos,oPiPphot,PtTh A*o, t'los,o|^ x) 

The full probability distribution can be derived from the phase-space distribution function 
/x which is the number of stars per unit velocity and per unit parallax, given by 

Up I) = fLip l)V{p) (24) 
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where V(p) is the volume per unit parallax at p (oc in three-dimensional space and oc p~^ 
for our zero-thickness disk), and is the number of stars per unit volume of phase space 
predicted by the model with parameters x: 

/x(P) ^tan, ^^los, = 1)U^{Vr, V^). (25) 

In the above equation /) is the number density of stars per unit area at (p, /) given 
by Equation ([2]) and ?7x is the number of stars per unit velocity predicted from the model 
parameters x at position {p, I) given by 

Uy,{vR,v^) = g{vR,VR^^,aR^^)g{v^,v^^^,a^^^). (26) 

Here g{y,y, cr) denotes the value at ?/ of a Gaussian distribution with mean y and dispersion 
0"^ and the quantities (I^_r,x, ^0,x) and {crR^y_, c"|,x) the mean velocity and velocity dispersion 
at parallax p and longitude / from the model with parameters x, given by Equations ([3]), 
(jlj), ([9]), and (|TOl) respectively. Finally, (ftan = A^/p, "yios) can be transformed to (f_R, f 0) for 
given (p, /), fview and Rq. 

In our experiment, we first observe Pphot for a random disk sample, with an error distri- 
bution e(pphotb) = g{Pphot,P,Sp,phot) about p. The distribution in (p,Pphot) of stars is given 
by 

S(p,Oe(Pphotb)^(p). (27) 

A subset of these stars is selected for our sample, with specified distribution Ai'sampie(Pphot, 
along a given line of sight, which can be related to the selection function, S'(pphot, (i-^-; ^^e 
probability of including a star in the survey at (pphot, 0) by 

POO 

A^sample(Pphot,0 = 'S'(Pphot,0 / ^{pJ)^{Pphot\p)V {p)dp. (28) 

Jo 

The total number of stars in the sample N is given by summing the number of stars toward 

roo 

Nl= iVsample(Pphot,Oc^Pphot, (29) 

^0 

and along all adopted lines of sight. 

Hence the full distribution of properties of the sample will depend on its intrinsic dis- 
tribution in phase-space /x, filtered by ^(pphot^O and convolved with appropriate error dis- 
tributions for the remaining observables, e(ptri|p) = g{ptTi,P,Sp,tTi), e(/^o|At) = 5'(/^o, ^/i) and 
€(■^108,0 1 "^los) = 5'('yios,o, wios, ^i;ioJ (as Outlined in Section [2.1.31) . The probability of finding a 
star in the survey is given by 



ot; Ptri; /^o; ^los, 
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S{PphotJ) f f dfJ^ f dViosfAP,f^/P,'"los,l)(^{Ptvi\p)(^{Pphot\p)^{fJ'o\f^)^{V\os,o\'"\os) 

J —oo P J —oo 

POO 

/ ^ipJ)Vip)eiptri\p)eiPphot\p)Pp,i,,vi,SP^f^o,Vios,o)dp, (30) 

^0 



Ni Jo 

_ 'S'(pphot,0 

Ni Jo 

where 



/ / U,,{fl/p, Vios)e{lIo\fJ')^{Vlos,o\Vlos) — dvios- (31) 
J —oo J —oo P 



Substituting in Equation fl23l) gives 

-P(j5phot,Ptri, Ho, Vlos,o\l, X) 



P{Ptri,IJ'0,Vlos ,o|Pphot) ^) X 



P(Pphotl^x) 

^ /o°° S(P. 0^(p)g(Ptrib)e(Pphotb)^p,M.^'ios(P^ /^o^ ^ios,o)^P ^22) 

/o°° C^Ptri JZo I-oo /o°° ^(P' 0^(p)e(Ptrib)e(Pphotb)^'p,M,^lo.(P, /^o, Vlos,o)dp 

Note that this expression is independent of our sample selection function ^(pphot, 
our analysis method leaves us free to choose a sample with arbitrary properties without 
biasing the results. Also, in the limit of negligible errors in the photometric parallax (i.e., 
e(Pphotb) = ^(Pphot -p)), Equation ([32]) simplifies to: 

-f (Ptri) A'-O) "^^loSjO 

|Pphot,^x) = e(ptribphot)-Pp,Ai,i>io,(Pphot ) (33) 

and our approach becomes insensitive to the underlying disk surface density distribution. 



2.2.3. Parameterizing the OUTPUT model 

The model distribution function /x (see Equation [2l]) is fully specified by the spatial 
number density, mean velocities and velocity dispersions of stars as a function of position in 
the disk. 

Our OUTPUT model can represent observations of axisymmetric motions using a total 
of 12 free parameters (11 free parameters when we fix Ro) to describe both /x and the 
transformation from physical to observed coordinates: (1) h — the scale length of the Galactic 
disk given in Equation ([5]); (2) (Tr^q and her (assuming the functional form for the radial 
velocity dispersion given in Eq. [H]); (3) the masses Mj {i = 1, 2, 8) within Galactic radius, 
Ri (from which the mass and its derivative DM /OR at any radius are found using cubic 
spline interpolation); and (4) Rq. We also analyze the sample assuming a known value 
Ro = 8 kpc in the following sections because we anticipate that it will be well-constrained by 
other observations (e.g., adaptive optics observations of stars around the black hole at the 
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Galactic center; lEisenhauer et al.ll2003l ). The effect of allowing Rq to be a free parameter or 



fixed is shown in Section 13.21 and 

All other quantities needed are derived from these parameters using the expressions 
given in Section r2.1.1[ 

Describing the spiral arms assuming m = 2 requires an additional four parameters: the 
constant c which is related to the radial wave number, k = c/R; the pattern speed flp] the 
arm's phase 0o and the amplitude of the spiral potential, The perturbations to the mean 
velocities and number density are calculated using Equations f fT^ . f|T5|) . and fl20|) . 

The 16 free parameters are listed in Table [H along with the INPUT values from which 
we derived our observed sample. 



3. Results 

In this section, we explore the accuracy of results recovered by applying the MCMC 
method to simulated observations of our disk model, with INPUT parameters given in Table 
[H For our standard sample, we look at disk M-giant stars, observed with fixed astrometric 
accuracies, as an example of a plausible near-future experiment that might be performed 
by a mission such as SIM Lite (see Section [3?T|) . M-giants are evolved, metal-rich stars and 
therefore typically relatively young (several Gyrs in age); this makes them good dynamical 
tracers of the mean disk potential. We then go on to examine how our results depend on 
sample size (Section [3. 2p . trigonometric and photometric accuracy (Sections 13. 31 and l 3.4p and 
disk coverage (Section l3.5p . 



3.1. The standard sample 



Our photometric sample is assumed to be composed of disk M-giants selected using the 
infra-red color in the 2MASS catalog. We assume that the intrinsic scatter in the absolute 
magnitudes of the M-giants around a mean of My = — 2 would result in a p hotometric 
parallax error of 15%, i.e., 5p,phot = O.lbp (as estimat ed bvlMaiewski et al.ll2003l ). Note that 
this scatter is in part due to metallicity differences (IChou et al.l 120071 ). and measuring the 
metallicity to about 0.3 dex would allow a parallax accuracy as good as ~10%. 

We select ~850 stars from our simulated M-giant survey to follow the distribution 
outlined in Section 12.1.31 and "observe" them with a trigonometric parallax accuracy of 



5. 



p,tri 



10 yuas. With these parameters, the point at which photometric rather than trigono- 
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metric parallaxes became more accurate would be at p ~ 67 /ias (~ 15 kpc). The proper 
motion accuracy is expected to scale as 

5^ = [0.235 + 0.634 x (5p,tri(/ias)] /lasyr"^ (34) 

from the SIM Global Astrometry Time Estimatoi@. We assume a constant error in the line-of- 
sight velocity 5„,ios=l kms~^ can be achieved from ground-based spectroscopic observations. 

Figure [1] illustrates the results of applying the MCMC recovery routine to this standard 
sample, with the analytical estimates constructed from INPUT and recovered OUTPUT 
parameters (listed in Table [1]) shown as dashed and solid lines, respectively. The best fit 
values at interpolated points are plotted as filled circles, with 1-a error bars estimated directly 
from the distribution of parameters in the MCMC. The figure indicates that, with this level 
of accuracy, we can recover the mass distribution between 4 and 20 kpc to within ~ 2% 
using a sample size < 10^. 

Figure [3] shows the full likelihood distribution of parameters. There are strong correla- 
tions between (a) Mj's and i^o^ due to the relation M = vl^^^RG'^; (b) Mj's at large R\ all 
Mi correlate with Rq simultaneously; (c) c and 0o: they define the phase of the spiral arms 
by Equations (fT^ and (fTSjl : (d) $a and VLp-. the amplitude of the mean velocity in the spiral 
arms is related to them by Equations (fT^ and (fT5l) . As noted in Section [2.2. H the steps 
in the MCMC were optimized using a preliminary run to take account of these correlations 
prior to our production runs. However, these contours in correlated parameters become 
longer and banana-shaped for small A^, large observational errors and narrow coverage of 
samples in space (e.g., small 0max defined in Section [3751) . and under these conditions the 
MCMC can fail to converge. 

In order to check for systematic biases in our methods as well as confirm the size of our 
error estimates. Figure [H repeats the middle left panel of Figured] for 10 runs of the MCMC 
(gray dots) applied to 10 independent samples of simulated stars. The sizes of the errors are 
not shown, but are similar to those shown in Figure [H Open circles and error bars indicate 
the mean AMj of the 10 runs and its estimated la error, and the solid line is the spline 
interpolation of the mean AMj. Only two of the means lie (a little) more than la away 
from the INPUT model (AMj = 0) and thus the scatter is consistent with the estimated 
statistical error of the means estimated from the 10 runs. In addition, the standard deviation 
of the OUTPUT mass parameters from their INPUT values, sdev = (0.20 ±0.02) x lO^^M© 
(calculated for the 10 runs over all eight points) is consistent with the mean of the errors 
estimated from the MCMC 5a/ = (0.21 ± 0.004) x IO^^Mq. Overall, these comparisons 



^http : //mscws4 . ipac . caltech. edu/simtools/portal/login/normal/l? 
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validate the error estimates derived from the MCMC (i.e., sdev ~ 5m) as well as the success 
of our method in modeling the adopted observing strategy. 

Note that the results in this section also confirm that this method does not require 
knowledge of the true spatial distribution of our sample stars — the algorithm we chose to 
pick the sample does not enter into the analysis. Indeed, the gray points in Figure [2] show 
that the intrinsic scatter in M-giant absolute magnitudes means that some of our sample 
lies outside of our intended survey region, i? = i^g = 20 kpc, where the last of our mass 
parameters is defined. A small systematic bias is apparent if we relax our requirement that 
a star's position defined by its photometric parallax lies well within this outer radius limit 
(recall, -Rmax = mm(18kpc,i?8 — O.lTrf)) and instead include stars whose observed Pphot 
places them all the way out to i? = 20 kpc. Our choice of keeping our samples to within 
= 18 kpc keeps this bias negligible. 

In each of the following subsections (and Figures [5] - [8]) , we repeat the comparison of 
sdev and 5m to check for systematic biases that may become apparent for samples observed 
under different conditions. 



Figure [5] plots sdev and 5 ]]/[ cLS cL function of the number of stars observed N, with 
all other properties of the sample maintained at their standard values (see Table [2] for full 
listing of errors on all parameters). The agreement of the black filled and gray open symbols 
demonstrates both the lack of systematic biases in our recovery algorithm and the success 
of the MCMC error estimates. The solid lines, representing the power law 



for fixing (fitting) Rq, confirm the expected A^~^/^ scaling of errors. The uncertainties in 
parameters for the case of fitting i?o are increased only by order unity compared to the case 
of fixing Rq. The uncertainty in Rq as a function of is given by the formula 



3.2. Dependence on 




(35) 




(36) 
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3.3. Dependence on Sp^tri 

Figure [6] and Table [3] summarize the results of repeating the analysis of Section 13.11 for 
MCMC runs based on the standard sample but with varying trigonometric parallax errors 
(5p tri (and related proper- motion accuracies — see Equation . 

The figure indicates that trends of the uncertainty Sm with 6p^tri can be roughly split 
into three regimes. As might be expected, 6m increases with dp^tri for 5p,tri less than ~ 10 
/ias and ^p^tri greater than ~ 200 /ias. However, 6m is almost constant for Sp^tri in the range 
10 — 200 /ias. This behavior can be understood by considering the importance of the sources 
of observational error in each of these three regimes. 

For 5p.tri < 10 yuas the trigonometric parallax is more accurate than the photometric parallax 
for most of the stars in the sample (those with p < 67 /ias and distances less than 15 
kpc), while the proper- motion (~ 6 /iasyr"^, corresponding to ftan < 1 kms~^ at < 20 
kpc) and line-of-sight velocity accuracies (1 kms~^) are much smaller than the scales 
of the velocity dispersions (~ 25 kms~^) of the population that they are trying to 
measure. Hence, the uncertainty scales with 5p tri 

For 10 /ias < Sp^tn < 200 /ias, the photometric parallax — held constant at 15% — provides 
stronger constraints on the results than the trigonometric parallax, and the proper 
motion and velocity accuracies are still too small to increase 6m- 

For 5p tri > 200 /ias, while the photometric parallax still provides 15% constraints on the 
distances to stars, the proper-motion error is increasing with ^p^tri beyond 0.21 masyr"^, 
or ~ 10 kms~^ for stars at 10 kpc. Hence, the accuracy with which the motions of 
distant stars can be determined is of the same order as their velocity dispersion and 
this now limits the accuracy with which the mass can be measured. 

We will discuss the implications of these trends for future surveys in Section |H 

The overlap of the black and gray points in Figure M once again confirms the lack of 
systematic biases in our recovery algorithm. 

3.4. Dependence on 5pphot 

We assume 15% photometric parallax accuracy in most of this work, i.e., 5p,phot = 
0.15p. We show the dependence of the accuracy of parameter estimates on various parallax 
accuracies, 5p,phot = O.lp, 0.15p and 0.2p for 5p tri = 10 and 1000 /tas in Figure [7] and Table HI 
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In general, we find that the mean error on parameters decreases shghtly as 5p,phot increases 
for 5p tri = 10 yuas, while the mean error is independent of 5p,phot for Sp,tTi = 1000 /ias. 

This counter-intuitive result — clearly in contradiction with the expectation that the er- 
rors on estimates should decrease as measurements become more accurate — can be attributed 
to the different nature of the samples in each case. Since the stars are selected using their 
photometrically estimated distances, (iphot (see Section [2.1.3p . the spatial distribution of the 
true positions of stars changes with ^p^phot^ the larger 5p,phot, the more the distribution of 
true distances is determined by the parameters of the disk rather than the parameters of the 
survey. 

For samples that are then observed with 5p tri = 10 /ias, the trigonometric parallax is 
more accurate than the photometric parallax for the majority of the stars {d < 15 kpc). 
Hence, even if 5p,phot increases, the accuracy of the dominant distance estimates in the 
analysis remains the same, while the area of the disk explored by the selected stars goes up 
and the net effect is an improvement in the errors on the parameters. 

For the same samples observed with 5p tri = 1000 //as, the photometric parallax is more 
accurate than the trigonometric parallax for the majority of the stars (> lOOpc). Now both 
the errors in the distance estimates and disk coverage increase with (5p,phot, and these effects 
exert competing influences on the parameter estimates. Hence the size of uncertainties on 
the parameters is largely independent of 5p phot- 

The only exception to these trends in errors is for the disk scale length h, where the 
mean error decreases as Sp^p^ot increases for both values of 5p,tri- This can be explained by 
recalling that the number density distribution only contributes to the likelihood function if 
the scattering of stars due to the uncertainty of Pphot is sufficiently large to sense the shape 
of the number density, by which the likelihood function is weighted — in the idealized case 
of zero errors, the method cannot constrain h at all (see Section 12.2.21 and Equations [23] - 
[30l) . This effect compounds rather than competes with the trend due to changing sample 
distributions as 5p,phot changes for both values of 6p^tri- 

To check how the results would be affected by systematic errors, we ran our analysis 
on samples constructed assuming photometric parallaxes that were 10% smaller or larger 
than the true parallax in addition to the random errors. The resultant Mj are systematically 
overestimated and underestimated by up to 10%, respectively. The parameters $a and h are 
also biased by 30%-50%. 
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3.5. Dependence on disk coverage and knowledge of Rq 

To check how the spatial distribution of our sample stars affects our recovery, we ran the 
MCMC for samples chosen with various values of 0max, the maximum absolute azimuthal 
angle of stars in the sample around the GC, i.e., |0| < (/'max (see Figure [2]). Here the angles 
(f) of sample stars are estimated using / and Pphot- In ah previous analyses we have not 
restricted (except that 

jnax ~ 133° in practice in our sample selection as seen in Figure 
[2]). The results are shown in Table [5] (with sample sizes = 850 and = 2000) and 
summarized in Figure [8] (for = 2000). As might be anticipated, the smaller the disk 
coverage (smaller 0max), the bigger the uncertainty of parameters — it is harder to be certain 
of the nonaxisymmetric features in the disk without a global view. For our particular disk 
model, it was necessary to cover more than 0max = 60° to recover parameters effectively. If 
0max was Sufficiently small, the MCMC failed to converge altogether (e.g., if 0max < 20° and 
< 10° with A^ = 850 and 2000, respectively, when Rq is fitted). 

Figures O and [S] and Tables [2] and O illustrate the results of experiments both for the case 
of fitting the distance Rq to the GC and the case of fixin g -Rn = 8 kpc assuming an accurate 



assessment of Rq has been made from other sources (e.g.. lEisenhauer et al.ll2003l ). If the disk 
is surveyed globahy (i.e., for large 0max), Ro can be recovered with a few percent accuracy 
(see Table [2] and E]) using the adopted samples, and the uncertainties in other parameters are 
increased only by order unity compared to the case of fixing Rq. However, for 0max < 60° 
the uncertainty in Rq and all other parameters increases even more dramatically with 0max 
compared to the examples where Rq was fixed. 



4. Discussion: implications for near-future surveys 
4.1. Astrometric surveys 

4-1.1. NASA's Space Interferometry Mission (SIM) Lite 

SIM Lite is a planned astrometric satellite using an optical interferometer, which will 
yield parallax errors as small as 4 /xas in wide angle mode. SIM Lite observes a parallax 
and a proper motion by repeatedly pointing at a target star and integrating as long as is 
necessary to get the required accuracy. For example. Figure [9] shows the expected mission 
time for stars requiring accuracies Sp^tri = 8, 10, 20 and 100 /ias observed in wide angle mode 
up to the limiting magnitude V = 20. [f| 



^http : //mscws4 . ipac . caltech.edu/ simtools/portal/login/nonnal/l? 
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Observers are allocated a specific amount of time to design their own experiments. 
Hence, for any disk study it is vital to find the optimal number, spatial coverage and re- 
quired trigonometric accuracies of target stars to recover the Galactic mass distribution most 
effectively. Figure [H] and Table [3] suggest that the uncertainties in our mass estimates are only 
weakly dependent on 6p^tn in the 10-100 /las regime. On the other hand. Figure [5] and Table 
|2] suggests that the uncertainty scales as \/l/N (as given in Equations [35] and [36]) . These 
results imply that our best strategy for studies of the mass distribution in the Galactic disk 
using SIM Lite (assuming that the astrometric errors are indeed statistical, not systematic) 
is to choose as large a sample as possible within the given observing time rather than using 
the time to get the best astrometric accuracies for a smaller number of targets. 

In principle, the number of stars observed by SIM Lite could be maximized by looking 
only at objects brighter than some limiting magnitude, V\. In order to make a realistic 
assessment of the extent to which such a magnitude limit would compromise our spatial 
coverage we need to also account for extinction in the disk plane. G. Zasowski, private 
communication, found typical V^-band extinctions in the Galactic disk to be: Ay ~ (1.5 — 
8) X rfkpc-i at |/| < 15°, Av ~ 0.8 x rfkpc"^ at |/| = 60° and Ay ~ 0.68 x rfkpc"^ at 
|/| = 90° and 180°. With our small sample size (A^ ~ 1000), we assume we can restrict 
our attention to stars that can be observed in low-extinction windows. Hence we adopt the 
lowest value. Ay ~ 1.5 x dkpc"^ for |/| < 15°. For |Z| > 15°, the typical values are linearly 
interpolated in /, and a fixed fraction Av,win/ ^v,typ of the estimated extinction is used. Figure 
[TU] illustrates the spatial coverage attainable for M-giants brighter than Vi = 16/18/20 and 
with Av^^in/ Av,typ = 0.5 (left-hand panel) and 0.8 (right-hand panel). While the extinction 
keeps us from seeing entirely across the Galactic disk in all cases, the requirement that 
0max > 60° in order to account for disk asymmetries (see Section ]33]) is still met, at least at 
small radii. 

Table [S] shows the number of M-giant stars that could be observed by SIM Lite within 
an allocation of 240 hr (motivated by the typical sizes of Key Project proposals for SIM 
PlanetQuest, the predecessor mission to the SIM Lite mission, which is similar, but based on 
a modified instrument architecture) for various Vy and Av^win/ Ay^typ, adopting the strategy 
outlined in previous sections of sampling equal numbers of stars in the 7 radial bins between 
the Ri in our model (see details in Section [2. 1.3p . The blank entries in the table correspond 
to cases where one or more of the radial bins contained no stars. Indeed, while the brighter 
limiting magnitudes do allow more stars to be observed (i.e., larger A^), the outer radial bins 
in these cases either contain no stars, or are only populated over a small range in 0. 

Overall, we find that the optimal sample has ~ 1000, Vi = 18 and 5p,tri = lOO/xas. 
Table [7] illustrates this by showing the uncertainty 6m for samples with 5p,tri = lOO/xas and 
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various V\. Even for the worst case of v4y,win/^y,typ = 1-0, the masses can be constrained 
with an uncertainty of 6m = 0.31 x IO^^Mq, which corresponds about 2% — 4% accuracy at 
R = 6- 20kpc. 

4.1.2. ESA 's GAIA satellite 

Unhke SIM Lite, GAIA is an all-sky survey satellite, where each star brighter than 
\^ = 20 is observed for the same amount of time. Concentrating again on the M-giants, 
this limiting magnitude means that the GAIA catalog will contain the stars in the disk 
as shown in Figure [101 although they will not be distributed uniformly in Galactocentric 
radius (because of the disk's intrinsic density gradient) and the variation of extinction along 
different lines of sight means that the depth of the sample will not be constant. In addition, 
the fainter stars will have less accurate astrometric measurements than the brighter ones 
(5ptri = 21 /ias for V = lb and 275 /xas for V = 20, with the corresponding proper-motion 
accuracies of 11 /lasyr"^ and 145 /xasyr^^ ). 

Figure [6] suggests that the target proper motion accuracies (better than ~ 100 /xasyr^^) 
are at the appropriate level to accurately assess the kinematical properties of disk stars when 
coupled with a photometric distance estimate. Nor do we expect the non-uniformity in the 
sample to introduce biases in results as our method is independent of the spatial distribution 
of the target stars. The gradient in the disk density coupled with the magnitude-limited 
nature of the survey means that the outer parts of the disk will be much more sparsely 
sampled and with larger error bars on the observations than the inner parts, so uncertainties 
in the mass estimates at large Galactocentric radii will increase correspondingly. In addition, 
since GAIA, like SIM, works in the optical, significant coverage beyond the Galactic center 
may be impossible due to extinction effects. However, these uncertainties can perhaps be 
offset by the sheer number of stars in the catalog. For example, the 2MASS catalog contains 
millions of M-giant candidates (i.e., in the color range 0.95 < J — K < 1.2) brighter than 
K = 14 (corresponding to V = 18) within 10 deg of the Galactic plane. While there will 
be significant contribution by bulge stars in this sample for Galactic longitudes |/| < 10°, 
they could be accounted for by adjusting the model to include this extra component. This 
argument assumes that the astrometric errors in GAIA are mainly statistical; given samples 
of this size, the biggest uncertainty with GAIA may be how large N can be before systematic 
errors begin to dominate the error budget. Combining the GAIA analysis with results from 
surveys that can fully cover the extent of the disk (e.g., SIM Lite; see Section 14.1.11) with 



http : //www ■ rssd . esa . int/index . php?proj ect=GAIA| ;page=Inf o_sheet s_overview 
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well-characterized systematics (e.g., VERA; see Section I4.1.3P should yield a global picture 
of the mass distribution within the Galactic disk. 



4.1.3. Radio VLBI arrays 



VERA jHonma et aPboOoh . VLBA jReidlboOsI : Iflachisuka et al.ll2009h . and EVN jRvd et al 



20081 ) are radio VLBI arrays that are conducting ~ lOyuas astrometric observations for water 
and/or methanol masers. Compared to our M-giant sample, water masers have the advan- 
tage as targets that most of them lie in star-forming regions very close to the Galactic plane. 
Hence they have rather lower velocity dispersion and fewer sources may be needed to ac- 
curately trace the rotation curve (and mass distribution). A possible disadvantage is that 
star-forming regions and therefore the maser targets are generated in spiral arms and thus 
tend to be found in a narrow range of azimuthal phases relative to the arms. In this case 
there may be larger uncertainties in and strong covariances between some of the parameters. 
Another disadvantage of the masers formed in young star-forming regions is that these tar- 
gets may not yet be dynamically homogenized to the disk, and so their motions may reflect 
the particular dynamics of their particular star-forming region, including such things as a 
vertex deviation, or other peculiar motions. 

The most exciting property of these samples is that they are already being produced. 
For example, 1 8 masers with pa rallax errors 6-80 /las by VLBA and VERA are already in 



the hterature (IReid et al.l l2009l ). and VERA is expected to observe approximately 1,000 
such sources over the next 10 years This sample will have observational error properties very 
similar to those assumed for our standard sample (Section 13. ip . except that photometric 
parallaxes are not available. Unlike our proposed M-giant survey, the distribution of the 
sources (mostly star-forming regions and some Mira variables) cannot be chosen arbitrarily, 
so the disk will not be uniformly sampled-in particular, water masers are more rare at large 
Galactocentric radii (> 10 kpc) and it may be hard to achieve full disk coverage. However, 
our results (e.g.. Figure E]) suggest that a sample of this size and level of accuracy could 
provide strong constraints on the mass distribution in the inner Galaxy. 

In addition to representing a significant step forward in measuring the Milky Way's 
rotation curve, these observations can serve as a vital cross check for future satellite surveys 
that rely on more sophisticated technology and hence may be more prone to unanticipated 
systematic biases. Furthermore, these radio VLBI observations do not suffer from dust 
extinction that the optical astrometry satellites GAIA and SIM Lite will. 
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4.2. Spectroscopic surveys 

High-resolution multi-object spectrographs are also currently under development. For 
example, the Apache Point Observatory Galactic Evolution Experiment (APOGEE — part 
of the Sloan Digital Sky Survey III project) will carry out a massive radial velocity survey 
starting in 2011. APOGEE will use a 300 fiber near-infrared (H-band) spectrograph with 
R = 22,500 — 25,000 to measure radial velocities to better than 0.5 kms^^ for more than 
10^ stars predominantly across the bulge and disk. The survey will take advantage of the 
low reddening in the near-infrared to reach stars throughout the Galactic disk and exploit 
the high resolution to make accurate estimates of spectroscopic parallaxes. Our results 
from Section 13.31 indicate that coupling these derived spectroscopic parallaxes with proper- 
motion measurements accurate to only ~ 1 masyr"^ could provide strong constraints on 
the mass distribution, especially given the large sample size. The Large Synoptic Survey 
Telescope could provide these proper motions for the third and fourth Galactic quadrants, 
down to a limiting magnitude of r = 24, allowing stars to be surveyed all the way across the 
Galactic disk. Partial coverage of the first and second quadrants can be achieved with similar 
accuracies (though brighter limiting magni tude) by corabinin g the Sloan Digital Sky Survey 



with the Palomar Observatory Sky Survey (IMunn et al.ll2004l ) . This approach would offer an 



alternative, independent assessment of the mass distribution to complement the astrometric 
measurements described above. 



5. Summary and Conclusions 

In this paper, we examined how accurately we might be able to recover the mass dis- 
tribution in the Galaxy (or more precisely, the gravitational force field in the Galactic disk) 
using current and near-future astrometric, photometric, and spectroscopic surveys of disk 
stars. We simulated observations of stars drawn from a simple model for the phase-space 
structure in an equilibrium, nonaxisymmetric disk and used a Markov Chain Monte Carlo 
approach to attempt to recover the model's parameters. Our formulation of this method 
relied on finding the parameter set of the model that maximizes the probability of stars 
in a survey having their observed trigonometric parallax, proper motion, and line-of-sight 
velocity given their measured photometric parallax. Hence, it is immune to biases in sample 
selection. Indeed, with a correct representation of the error distributions for each observable 
and of the non-axisymmetric components of the mass distribution, no systematic errors were 
evident in our approach for samples of 100s- 1000s of stars combining observed trigonometric 
parallaxes with errors in the range 1 /ias to 2 mas (corresponding to proper motion errors 
assumed to be in the range 0.9 /iasyr"^ to 1.2 masyr"^), photometric parallaxes known at 
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the 10%-20% level and line-of-sight velocity measurements accurate to 1 kms ^. 

The presence of non-axisymmetric features in the disk means that a precise mapping of 
the Galactic mass distribution will require a survey on global scales. Even in our simplified 
case with a single two- armed spiral pattern, restricting our survey to a limited range of 
Galactic longitude significantly reduced the accuracy of our results. 

However, given such a global disk survey, we found that we could recover the mass 
profile in the range 4-20 kpc with a few percent accuracies using a variety of approaches. 
If lias trigonometric parallaxes are available (with associated proper-motion measurements 
and 1 kms~^ line-of-sight velocities), then this accuracy is feasible with a survey as small as 
a few hundred stars. Once trigonometric parallax errors exceed 10 //as, the same accuracy 
can be achieved by supplementing the trigonometric parallaxes with photometric parallaxes 
accurate to 10%-20% and adopting a sample of thousands of stars, so long as proper-motion 
errors remain below a level of few hundred //asyr"^. If proper-motion errors are of order a 
few masyr"^, then larger samples are needed in compensation. 

We also found we could measure the mass distribution even in the absence of an accurate 
assessment of the distance to the Galactic center, i?o- Including Rq as a free parameter did 
increase our uncertainties by a factor of 2, but also allowed us to measure this distance 
with comparable accuracy to the mass distribution itself (i.e., a few percent for the samples 
discussed above). 

We conclude that, whether one or all of the future surveys (e.g., SIM Lite, GAIA, VERA 
and APOGEE) are completed, a significant step forward in our understanding of the Galactic 
mass distribution (i.e., an assessment of the force field in the Galactic disk at the 1% level) 
is on the horizon, as well as detailed insights into disk dynamics. 

These conclusions are based on the assumption that the deviations from a smooth 
axisymmetric model of the gravitational field in the disk can be modeled as a grand-design 
spiral pattern with a specified form (logarithmic spiral) and a well-defined pattern speed, 
without resonances within the disk. Further work is required to understand how relaxing 
these assumptions would affect the accuracy of astrometric disk surveys. Natural directions 
for future works are including vertical motions of sample stars as well as a bulge component, 
and more general models of spiral structure. 

We acknowledge R. Patterson and J. Carlin for useful assistance with calculations re- 
garding SIM and 2MASS. We thank L. S. Eycr for providing important insights into GAIA. 
We are grateful to M. Honma for the comments about VERA. We thank G. Zasowski for help 
with the dust extinction estimates used here. This work is supported by NASA/ JPL contract 
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Table 1. Parameters with INPUT values and an example OUTPUT. 



Parameter 


-Ri(kpc) 


INPUT 


OUTPUT 


1(7 Error 




4.00 


2.37 


2.24 


0.15 


M2(lOiOM0) 


6.29 


4.94 


4.90 


0.10 




8.57 


7.86 


7.77 


0.12 


M4(10i°Afo) 


10.86 


10.92 


10.76 


0.16 


M5(lOiOM0) 


13.14 


13.99 


13.83 


0.17 




15.43 


16.99 


17.10 


0.27 




17.71 


19.90 


19.72 


0.30 




20.00 


22.68 


22.77 


0.33 


(TR,0(kms-^) 




25.00 


24.64 


0.64 


00 (rad) 




1.83 


1.85 


0.07 


c (rad) 




7.46 


7.51 


0.17 


Qp (kms~^ kpc^^) 




1.40 


1.12 


0.37 


$a(kms"^)^ 




200.0 


208.3 


32.8 


K (kpc) 




3.00 


3.14 


0.10 


h (kpc) 




3.00 


3.36 


0.35 


Ro (kpc) 




8.00 


7.97 


0.12 



Note. — Recovered OUTPUT parameters and la errors for a 
survey with (5p,phot = 0.15p and Sp^tii = 10 /xas and N — 850 stars. 
Mi is the mass within Ri. 
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Table 2. Mean errors of parameters Versus number of stars N. 



N 


M 






c 


S7p 




her 


h 


Ro 


500 


0.17 


0.79 


0.08 


0.23 


0.54 


48.2 


0.12 


0.44 




850 


0.13 


0.60 


0.06 


0.17 


0.41 


36.7 


0.09 


0.29 




2000 


0.08 


0.39 


0.04 


0.10 


0.30 


27.5 


0.06 


0.19 




4000 


0.06 


0.28 


0.02 


0.07 


0.20 


19.1 


0.04 


0.12 




8000 


0.04 


0.20 


0.02 


0.05 


0.13 


12.4 


0.03 


0.09 




500 


0.27 


0.88 


0.09 


0.22 


0.52 


44.8 


0.12 


0.43 


0.16 


850 


0.21 


0.66 


0.07 


0.17 


0.42 


37.5 


0.09 


0.32 


0.12 


2000 


0.13 


0.42 


0.05 


0.11 


0.30 


27.2 


0.06 


0.19 


0.08 


4000 


0.09 


0.30 


0.03 


0.07 


0.21 


19.1 


0.04 


0.13 


0.05 


8000 


0.07 


0.21 


0.02 


0.05 


0.13 


11.7 


0.03 


0.09 


0.04 



Note. — The means of the estimated error of parameters are calculated 
from 10 runs of MCMC for the case of Sp^tri = 10 A^a-s, either fixing Rq = 8 
kpc (shown as — ) or fitting Rq. Units are same as in Table [T] For M, the 
mean is also taken over i =1-8. See the notes for Table [Tl 



Table 3. Mean errors of parameters with various trigonometric parallax errors 5p^tr\- 



5p,tTi 


M 




<t>(i 




c 






ha 


h 


1 


0.10 


0.58 


0.05 


0. 


,13 


0.33 


32.8 


0.07 


0.22 


2 


0.10 


0.57 


0.05 


0. 


,14 


0.34 


33.5 


0.08 


0.26 


4 


0.11 


0.58 


0.05 


0. 


14 


0.40 


37.3 


0.08 


0.28 


8 


0.12 


0.59 


0.06 


0. 


,16 


0.41 


37.3 


0.09 


0.29 


10 


0.13 


0.60 


0.06 


0. 


,17 


0.41 


36.7 


0.09 


0.29 


20 


0.14 


0.61 


0.06 


0. 


,18 


0.43 


38.2 


0.09 


0.32 


50 


0.15 


0.63 


0.06 


0. 


,18 


0.44 


38.2 


0.10 


0.35 


100 


0.15 


0.65 


0.06 


0. 


,18 


0.45 


39.0 


0.10 


0.37 


200 


0.16 


0.65 


0.06 


0. 


,19 


0.45 


38.7 


0.11 


0.38 


500 


0.18 


0.67 


0.06 


0. 


,20 


0.47 


40.3 


0.12 


0.40 


1000 


0.23 


0.68 


0.07 


0. 


,22 


0.50 


41.8 


0.14 


0.54 


2000 


0.31 


0.74 


0.07 


0. 


,23 


0.56 


46.1 


0.15 


0.57 



Note. — For A'^ = 850 and fixing Ro = 8 kpc. See the notes for 
Tabled 
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Table 4. Mean errors of parameters with various photometric parallax errors (5p,phot- 



<5p,tri 




M 






c 






ha 


h 


10 


0.10 


0.14 


0.60 


0.06 


0.17 


0.48 


42.0 


0.09 


0.45 


10 


0.15 


0.13 


0.60 


0.06 


0.17 


0.41 


36.7 


0.09 


0.29 


10 


0.20 


0.12 


0.59 


0.05 


0.15 


0.39 


37.6 


0.09 


0.23 


1000 


0.10 


0.23 


0.66 


0.06 


0.21 


0.51 


40.5 


0.12 


1.07 


1000 


0.15 


0.23 


0.68 


0.07 


0.22 


0.50 


41.8 


0.14 


0.54 


1000 


0.20 


0.23 


0.70 


0.06 


0.20 


0.47 


41.0 


0.14 


0.37 



Note. — For A'^ = 850 and fixing Ro = 8 kpc. (5p tri is in fias. See the notes 
for Table [2] 



Table 5. Mean errors of parameters with various angular coverages 0] 



N I 


^max 


M 




c 


4>o 




'S'a 


ha 


h 


Ro 


850 


133 


0.13 


0.60 


0.06 


0.17 


0.41 


36.7 


0.09 


0.29 




850 


90 


0.13 


0.61 


0.06 


0.17 


0.41 


37.6 


0.09 


0.29 




850 


60 


0.14 


0.58 


0.05 


0.17 


0.39 


37.4 


0.09 


0.36 




850 


30 


0.18 


0.59 


0.06 


0.20 


0.37 


35.6 


0.08 


0.48 




850 


20 


0.22 


0.58 


0.07 


0.22 


0.42 


39.8 


0.08 


0.46 




850 


10 


0.40 


0.58 


0.06 


0.20 


0.45 


43.4 


0.08 


0.57 




850 


133 


0.21 


0.66 


0.07 


0.17 


0.42 


37.5 


0.09 


0.32 


0.12 


850 


90 


0.25 


0.68 


0.08 


0.17 


0.41 


37.7 


0.09 


0.33 


0.16 


850 


60 


0.37 


0.72 


0.11 


0.19 


0.41 


39.6 


0.09 


0.38 


0.24 


850 


30 


0.79 


1.15 


0.24 


0.34 


0.44 


41.4 


0.08 


0.40 


0.56 


2000 


133 


0.08 


0.39 


0.04 


0.10 


0.30 


27.5 


0.06 


0.19 




2000 


90 


0.08 


0.39 


0.04 


0.11 


0.28 


26.7 


0.06 


0.21 




2000 


60 


0.09 


0.39 


0.04 


0.12 


0.28 


26.5 


0.06 


0.23 




2000 


30 


0.12 


0.39 


0.04 


0.13 


0.30 


29.5 


0.06 


0.25 




2000 


20 


0.15 


0.37 


0.04 


0.13 


0.29 


29.1 


0.05 


0.27 




2000 


10 


0.25 


0.38 


0.04 


0.14 


0.31 


29.6 


0.05 


0.29 




2000 


133 


0.13 


0.42 


0.05 


0.11 


0.30 


27.2 


0.06 


0.19 


0.08 


2000 


90 


0.16 


0.42 


0.05 


0.11 


0.26 


24.5 


0.06 


0.23 


0.10 


2000 


60 


0.23 


0.59 


0.07 


0.13 


0.28 


28.0 


0.06 


0.26 


0.16 


2000 


30 


0.66 


0.61 


0.22 


0.31 


0.34 


30.6 


0.05 


0.27 


0.49 


2000 


20 


0.92 


0.61 


0.29 


0.38 


0.36 


33.6 


0.05 


0.31 


0.67 



Note. — The case of <5p t^; = 10 /^as. MCMC fails to converge for 0max < 20° 
and < 10° for N = 850 and 2000, respectively, in the case where Ro is fitted. "— " 
means that _Ro is fixed. See the notes for Table |2] 
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Table 6. Number of observable stars by SIM Lite in 240 hr with various extinctions and 

limiting magnitude. 





^V,win 








Ar(yi[mag]) 








{fias) 




^1=20 


Vi=19 


Vi=18 


Vi=17 


Vi=16 


yi=15 


yi=i4 


o 
O 


i.U 


0/1 














o 
o 


U.o 


O / 


/ U 












Q 
O 


U.o 


vo 
/ Z 


1 on 


ZLo 


ooD 








Q 

o 


U.o 


iz / 




QOQ 

oZo 


/I c;o 
40o 


Doi 


OOO 




Q 

o 


n n 
U.U 


1 r / 


oOU 


ODD 


III 


yo / 


lUDy 


1 1 Q 1 


10 


1.0 


45 














10 


0.8 


69 


128 












10 


0.5 


132 


225 


354 


521 








10 


0.3 


224 


355 


505 


669 


853 


1031 




10 


0.0 


304 


550 


779 


969 


1108 


1182 


1215 


20 


1.0 


119 














20 


0.8 


172 


296 












20 


0.5 


305 


470 


656 


864 








20 


0.3 


474 


664 


835 


996 


1126 


1212 




20 


0.0 


584 


856 


1030 


1147 


1211 


1235 


1237 


100 


1.0 


295 














100 


0.8 


399 


594 












100 


0.5 


613 


807 


974 


1108 








100 


0.3 


808 


969 


1083 


1163 


1209 


1233 




100 


0.0 


900 


1079 


1166 


1212 


1231 


1237 


1237 



Note. — The numbers of observable stars A^(yi[mag]) are estimated by using 
the required mission time for SIM Lite given by Figure |9] A factor j4v,win/^v,typ 
indicates the fraction of V-band extinction at low extinction window relative to the 
typical extinction in that direction (see Section [4.1.1l l. Vi represents V-band limiting 
magnitude to select sample stars. " — " means that stars can not be sampled in all 
Ri bins with that VJ. 



Table 7. 6m by SIM Lite in 240 hr with various extinctions and limiting magnitude. 



<5p,tri 






<5m 




(^as) 


Av,typ 


Vi=20 


Vi=19 Vi=18 


Vi=17 


100 


1.0 


0.311 






100 


0.8 


0.238 


0.216 — 




100 


0.5 


0.178 


0.161 0.149 


0.156 



Note. — are estimated for the cases given in Table 
[6]with 5p^tri = 100 /xas and fixing Rq. The survey regions 
are as shown in Figure [TOl The numbers of sample stars 
are as in Table |6] 



5 10 15 20 5 10 15 20 



R (kpc) R (kpc) 




R (kpc) R (kpc) 

Fig. 1. — INPUT (dashed lines) and recovered OUTPUT (solid lines with ±lcr dispersion) 
models as a function of the Galactocentric radius R. The OUTPUT is derived for samples 
of = 850 stars, observed with assumed photometric parallax and trigonometric parallax 
accuracies of (5p,phot = 0.15p and 5p,tri = 10 /ias. The model curves are interpolated from 
the mass parameters, Mj (indicated for the OUTPUT as points with la error bars). The 
oscillations in vji are due to spiral arms. Top-left: mass within R, M{R). Middle left: 
residual from the INPUT mass, AM{R). Top right: azimuthal velocity f^. Middle right: 
residual from the INPUT circular velocity, Afcirc- Bottom left: radial velocity dispersion an, 
which is estimated from Equation using OUTPUT ctr,© and her- Bottom right: radial 
velocity vr. The values of f,^ and vr for stars with = ±7r/16 are plotted as light gray dots. 
In some plots the error bars are too small to be resolved and the dashed and solid lines lie 
on top of one another. We can recover the mass to within ~2% for i?=4-20 kpc. 
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X (kpc) 

Fig. 2. — Spatial distribution of selected stars on the Galactic plane seen from the North 
Galactic Pole. The large gray dots indicate their intrinsic positions and the small black dots 
represent their positions derived from their "observed" photometric parallax. The sample is 
chosen from their observed positions so that there are equal numbers of stars in the seven 
bins defined by the Galactocentric radii Ri (Section I2.1.3p . The symbols + and indicate 
the GC and the Sun. We use stars (small black dots) between 4 (inner thick circle) and 18 
kpc from the GC and < 20 kpc from the Sun (dashed circle). We also select only stars within 
20kpc-0.17(i from the GC (short-long dashed line) to reduce the number of stars that leak 
outside the outermost model radius Rg = 20 kpc due to errors in the photometric parallax 
(see details in Section [2.1.31) . Stars that lie behind the GC with |/| < 30° (dashed lines) 
are not used. The peak of the spiral arm potential is also drawn. The angle is defined as 
shown by the arrow. 
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Fig. 3. — Likelihood distribution of parameters for the MCMC run shown in Figure [Tj with 
= 850, 5p,phot = 0.15p and Sp^tri = 10 /las. The dots indicate the maximum hkehhood val- 
ues. The crosses indicate the INPUT values. The contours indicate 68% and 95% confidence 
intervals. The histograms represent the projected likelihood distribution for each parameter, 
where the vertical lines indicate INPUT values. 
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Fig. 4. — Residuals from INPUT mass, AMj, of 10 MCMC runs as a function of R (gray 
dots), as in the middle left panel of Figure ^ N = 850, Sp^p^ot = 0.15p and 5p,tri = 10 
/xas. The distance to the GC, Rq, is one of the fitting parameters. Here, each run uses an 
independent sample of simulated stars. The open circles with error bars indicate the mean 
of 10 runs with la errors. The solid line is the interpolation of the mean AMj. One can see 
that any systematic differences are within the error, ~ 0.1 x lO^^M©. 
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500 1000 2000 5000 10^ 

N 

Fig. 5. — Mean of errors in Mj calculated by the MCMC {Sm'- filled symbols) and the 
standard deviation of estimates from an ensemble of 10 runs (sdev: open symbols with error 
bars) as a function of the number of stars, N, for 5p,phot = 0.15p, Sp^tri — 10 //as. The squares 
denote runs in which Rq was fitted and the circles denote runs in which Rq was fixed at 8 
kpc. Here the mean oi Sm and sdev is taken over i = 1-8 and over 10 runs of MCMC. The 
solid fines indicate the scaling SM/{10^Mq x y/500/N) = 1.7 and 2.7, respectively. 
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Fig. 6. — Mean of estimated errors in Mj and the standard deviation from the true Mj 
as a function of trigonometric parallax error for sample size N — 850, 5p,phot = 0.15p and 
<^p,tri = 10 /xas. In these runs Rq is fixed at 8 kpc. 
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Fig. 7. — Mean of estimated errors in Mi (filled symbols) and the standard deviation from 
the true Mj (open symbols) as a function of relative photometric parallax error 5p,phot/p for 
= 850, (5p_tri = 10 /ias (circle) and 6p^tri = 1000 /ias (square). In these runs Rq is fixed at 
8 kpc. Note that in some cases the error in derived quantity goes down as the observational 
error goes up (see discussion in Section [3^ . 
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Fig. 8. — Mean estimated error in Mj (filled symbols) and the standard deviation from the 
true Mi (open symbols) for N — 2000 as a function of 0max- The circles and squares are for 
the case of fixing it!o = 8 kpc and fitting Rq, respectively. The MCMC fails to converge for 
0max < 10° when i?o is fitted. Covering ^max > 60° is required to recover models efficiently, 
especially for the case of fitting Rq. 
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Fig. 9. — SIM Lite's required mission time to measure the parallax of a star with 8, 10, 20 
and 100 //as accuracy as a function of V-band magnitude of the star. The top horizontal 
axis indicates the distance to the target with an absolute magnitude of My — —2 mag (e.g., 
M giants). 
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Fig. 10. — Same as Figure [21 but showing the maximum survey region with the V^-band 
hmiting magnitude VJ = 16, 18 and 20 mag (the thick gray hues) for the case that we observe 
M-giants with an absolute magnitude of My = —2 through the low extinction windows with 
^y,wm/^y,typ =0.5 (left panel) and 0.8 (right panel) (see details in Section 14.1.11) . 



